# JoVE Jidesheng Snake Pills meta-analysis
# Source-verified R re-analysis script
# Requires: metafor, readxl
library(metafor)
library(readxl)

data_file <- "Supplementary_Dataset.xlsx"
output_dir <- "R_outputs"

run_mkh <- function(yi, vi, slab, label) {
  fit_z <- rma.uni(yi=yi, vi=vi, method="REML", test="z", slab=slab)
  fit_hk <- rma.uni(yi=yi, vi=vi, method="REML", test="knha", slab=slab)
  k <- length(yi)
  se_mkh <- max(as.numeric(fit_z$se), as.numeric(fit_hk$se))
  crit <- qt(0.975, df=k-1)
  est <- as.numeric(fit_z$b)
  data.frame(outcome=label, k=k, estimate=est, se_mkh=se_mkh,
             ci_lb=est-crit*se_mkh, ci_ub=est+crit*se_mkh,
             tau2=fit_z$tau2, I2=fit_z$I2, QE=fit_z$QE, QEp=fit_z$QEp)
}

dir.create(output_dir, showWarnings=FALSE)
dir.create(file.path(output_dir,"figures"), recursive=TRUE, showWarnings=FALSE)

# 1) Response: source-reported effects only; no primary pooled estimate
dich <- read_excel(data_file, sheet="R_Data_Dichotomous")
src <- subset(dich, endpoint_version=="source_reported")
src$yi <- log((src$events_t/src$n_t)/(src$events_c/src$n_c))
src$vi <- 1/src$events_t - 1/src$n_t + 1/src$events_c - 1/src$n_c
src$RR <- exp(src$yi)
src$ci_lb <- exp(src$yi - 1.96*sqrt(src$vi))
src$ci_ub <- exp(src$yi + 1.96*sqrt(src$vi))
write.csv(src, file.path(output_dir,"response_source_specific.csv"), row.names=FALSE)

# Post hoc harmonized any-improvement sensitivity only
harm <- subset(dich, harmonized_sensitivity=="YES")
harm_es <- escalc(measure="RR", ai=events_t, n1i=n_t, ci=events_c, n2i=n_c, data=harm, append=TRUE)
harm_re <- run_mkh(harm_es$yi, harm_es$vi, harm_es$study, "Harmonized any-improvement RR")
harm_re$estimate <- exp(harm_re$estimate)
harm_re$ci_lb <- exp(harm_re$ci_lb)
harm_re$ci_ub <- exp(harm_re$ci_ub)
write.csv(harm_re, file.path(output_dir,"response_harmonized_REML_mKH.csv"), row.names=FALSE)
harm_fe <- rma.uni(yi=harm_es$yi, vi=harm_es$vi, method="FE", slab=harm_es$study)
write.csv(data.frame(estimate=exp(as.numeric(harm_fe$b)),ci_lb=exp(harm_fe$ci.lb),ci_ub=exp(harm_fe$ci.ub),
                     I2=harm_fe$I2,QE=harm_fe$QE,QEp=harm_fe$QEp),
          file.path(output_dir,"response_harmonized_fixed_sensitivity.csv"), row.names=FALSE)

# 2) Continuous outcomes
cont <- read_excel(data_file, sheet="R_Data_Continuous")

sw <- subset(cont, outcome=="swelling_resolution_days" & primary_pool=="YES")
sw_es <- escalc(measure="MD",m1i=mean_t,sd1i=sd_t,n1i=n_t,m2i=mean_c,sd2i=sd_c,n2i=n_c,data=sw,append=TRUE)
sw_re <- run_mkh(sw_es$yi, sw_es$vi, sw_es$study, "Swelling-resolution time (days)")
write.csv(sw_re,file.path(output_dir,"swelling_REML_mKH.csv"),row.names=FALSE)
sw_fe <- rma.uni(yi=sw_es$yi,vi=sw_es$vi,method="FE",slab=sw_es$study)
write.csv(data.frame(estimate=as.numeric(sw_fe$b),ci_lb=sw_fe$ci.lb,ci_ub=sw_fe$ci.ub,I2=sw_fe$I2,QE=sw_fe$QE,QEp=sw_fe$QEp),
          file.path(output_dir,"swelling_fixed_sensitivity.csv"),row.names=FALSE)

loo <- do.call(rbind,lapply(seq_len(nrow(sw_es)),function(i){
  tmp <- sw_es[-i,]
  z <- run_mkh(tmp$yi,tmp$vi,tmp$study,paste0("Exclude ",sw_es$study[i]))
  z$excluded <- sw_es$study[i]
  z
}))
write.csv(loo,file.path(output_dir,"swelling_leave_one_out_REML_mKH.csv"),row.names=FALSE)

sw_no_li <- subset(sw_es,study!="Li BL 2022")
write.csv(run_mkh(sw_no_li$yi,sw_no_li$vi,sw_no_li$study,"Swelling excluding Li BL 2022"),
          file.path(output_dir,"swelling_exclude_Li_population_sensitivity.csv"),row.names=FALSE)

hs <- subset(cont,outcome=="hospital_stay_days" & primary_pool=="YES")
hs_es <- escalc(measure="MD",m1i=mean_t,sd1i=sd_t,n1i=n_t,m2i=mean_c,sd2i=sd_c,n2i=n_c,data=hs,append=TRUE)
write.csv(run_mkh(hs_es$yi,hs_es$vi,hs_es$study,"Hospital stay (days)"),
          file.path(output_dir,"hospital_REML_mKH.csv"),row.names=FALSE)
hs_fe <- rma.uni(yi=hs_es$yi,vi=hs_es$vi,method="FE",slab=hs_es$study)
write.csv(data.frame(estimate=as.numeric(hs_fe$b),ci_lb=hs_fe$ci.lb,ci_ub=hs_fe$ci.ub,I2=hs_fe$I2,QE=hs_fe$QE,QEp=hs_fe$QEp),
          file.path(output_dir,"hospital_fixed_sensitivity.csv"),row.names=FALSE)

# 3) CRP: narrative only — no primary raw-MD meta-analysis
crp <- subset(cont,outcome=="crp_post")
write.csv(crp,file.path(output_dir,"CRP_narrative_source_values.csv"),row.names=FALSE)

writeLines("Formal publication-bias testing was not performed because fewer than 10 studies contributed to each pooled outcome; publication bias was considered not assessable.",
           file.path(output_dir,"publication_bias_statement.txt"))

cat("Source-verified re-analysis complete.\n")

# Running this file reproduces the analysis and all five figures from the files
# supplied in the same folder. Set JOVE_ANALYSIS_ONLY <- TRUE before sourcing
# this script when analysis outputs alone are required.
if (!exists("JOVE_ANALYSIS_ONLY") || !isTRUE(JOVE_ANALYSIS_ONLY)) {
  source("Figure_1_Source.R")
  source("Figures_2_to_5_Source.R")
}